library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(gtable)
library(gtsummary)
## Warning: package 'gtsummary' was built under R version 4.1.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.1.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(knitr)
library(car)
## Warning: package 'car' was built under R version 4.1.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(boot)
## Warning: package 'boot' was built under R version 4.1.3
## 
## Attaching package: 'boot'
## The following object is masked from 'package:car':
## 
##     logit
data <- readRDS("Dataset.RData")
data$y = 0
data$y[data$survival == 'Long'] = 1
# immune selected genes MC1R, AKT1, NAMPT, TNFRSF19, PAK2, SOS1, B2M 
l1 <- glm(data$y ~ data$MC1R + as.factor(data$stage), family = 'binomial')
l1$coefficients[2]
##  data$MC1R 
## -0.9632995
l2 <- glm(data$y ~ data$AKT1 + as.factor(data$stage), family = 'binomial')
l2$coefficients[2]
## data$AKT1 
## -0.917543
l3 <- glm(data$y ~ data$NAMPT + as.factor(data$stage), family = 'binomial')
l3$coefficients[2]
## data$NAMPT 
##  0.1939064
l4 <- glm(data$y ~ data$TNFRSF19 + as.factor(data$stage), family = 'binomial')
l4$coefficients[2]
## data$TNFRSF19 
##       0.23159
l5 <- glm(data$y ~ data$PAK2 + as.factor(data$stage), family = 'binomial')
l5$coefficients[2]
## data$PAK2 
## 0.5984475
l6 <- glm(data$y ~ data$SOS1 + as.factor(data$stage), family = 'binomial')
l6$coefficients[2]
## data$SOS1 
## 0.2604371
l7 <- glm(data$y ~ data$B2M + as.factor(data$stage), family = 'binomial')
l7$coefficients[2]
##  data$B2M 
## 0.3445805
dd = data[c('MC1R', 'AKT1', 'NAMPT', 'TNFRSF19', 'PAK2', 'SOS1', 'B2M')] 
library(plotly)
## Warning: package 'plotly' was built under R version 4.1.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.1.3
library(heatmaply)
## Warning: package 'heatmaply' was built under R version 4.1.3
## Loading required package: viridis
## Loading required package: viridisLite
## 
## ======================
## Welcome to heatmaply version 1.3.0
## 
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
## 
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## You may ask questions at stackoverflow, use the r and heatmaply tags: 
##   https://stackoverflow.com/questions/tagged/heatmaply
## ======================
heatmaply_cor(cor(dd))
cor(dd)
##                MC1R       AKT1      NAMPT    TNFRSF19       PAK2        SOS1
## MC1R      1.0000000  0.6166103 -0.1633626 -0.13110913 -0.3211794 -0.35683222
## AKT1      0.6166103  1.0000000 -0.2861000 -0.24337733 -0.4490884 -0.37407132
## NAMPT    -0.1633626 -0.2861000  1.0000000  0.23196342  0.1865976  0.11443004
## TNFRSF19 -0.1311091 -0.2433773  0.2319634  1.00000000  0.2679136  0.04096849
## PAK2     -0.3211794 -0.4490884  0.1865976  0.26791357  1.0000000  0.21673064
## SOS1     -0.3568322 -0.3740713  0.1144300  0.04096849  0.2167306  1.00000000
## B2M      -0.2259642 -0.1928459  0.2008789 -0.08569906  0.1378465  0.16882559
##                  B2M
## MC1R     -0.22596423
## AKT1     -0.19284591
## NAMPT     0.20087893
## TNFRSF19 -0.08569906
## PAK2      0.13784648
## SOS1      0.16882559
## B2M       1.00000000
dd1 = as.data.frame(dd)
dd1$y = data$y
dd1$age_c = data$age - mean(data$age)
dd1$sex = 0
dd1$sex[data$gender == 'female'] = 1
dd1$pathStage = as.factor(data$stage)
model = glm(y ~ MC1R + AKT1 + NAMPT + TNFRSF19 + PAK2 + SOS1 + B2M + pathStage, data = dd1, family = 'binomial')
vif(model)
##               GVIF Df GVIF^(1/(2*Df))
## MC1R      1.199758  1        1.095335
## AKT1      1.435392  1        1.198078
## NAMPT     1.127966  1        1.062058
## TNFRSF19  1.183216  1        1.087757
## PAK2      1.197908  1        1.094490
## SOS1      1.192218  1        1.091888
## B2M       1.131922  1        1.063918
## pathStage 1.134084  3        1.021192
# selected top genes: MC1R, AKT1, PAK2, SOS1, B2M
# model 1: basline + selected top genes

model1 <- glm(y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + pathStage, data = dd1, family = 'binomial')
summary(model1)
## 
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + 
##     pathStage, family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5591  -0.9368   0.2876   0.9012   2.0907  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.83506    0.46974   3.907 9.36e-05 ***
## MC1R         -0.51267    0.30252  -1.695  0.09014 .  
## AKT1         -0.60824    0.24547  -2.478  0.01322 *  
## PAK2          0.32723    0.18588   1.760  0.07834 .  
## SOS1         -0.10986    0.17285  -0.636  0.52506    
## B2M           0.22297    0.15674   1.423  0.15488    
## sex          -0.23477    0.31377  -0.748  0.45432    
## age_c        -0.02149    0.01077  -1.995  0.04600 *  
## pathStage2   -1.64327    0.52583  -3.125  0.00178 ** 
## pathStage3   -2.33998    0.50959  -4.592 4.39e-06 ***
## pathStage4  -18.48491  882.74369  -0.021  0.98329    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 331.89  on 239  degrees of freedom
## Residual deviance: 251.47  on 229  degrees of freedom
## AIC: 273.47
## 
## Number of Fisher Scoring iterations: 13
model1.diag <- glm.diag(model1)
data.frame("cooks" = model1.diag$cook, "leverage" = model1.diag$h) %>%
  pivot_longer(everything(), names_to = "measure") %>%
  mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
  ggplot(aes(x = seq(0.5,240,by=0.5))) +
  geom_line(aes(y = value, col = measure)) +
  geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
  theme_classic() +
  facet_wrap(~measure, scales = "free_y", nrow = 2) +
  scale_color_manual(values = c("orange","cornflowerblue")) +
  labs(x = "observation")

p_res1 <- residuals(model1, type = "pearson")
d_res1 <- residuals(model1, type = "deviance")
data.frame("pearson" = p_res1, "deviance" = d_res1) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation")

data.frame("pearson" = model1.diag$rp, "deviance" = model1.diag$rd) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
  labs(x = "observation", y = "standardized values")

#kable(data.frame("pearson" = p_res1, "deviance" = d_res1,"leverage" = model1.diag$h, "cooks" = model1.diag$cook,"p_standard" = model1.diag$rp, "d_standard" = model1.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
# model 2: baseline + selected top genes + immune resposne genes
dd1$CD86 = data$CD86
dd1$DDX58 = data$DDX58
dd1$KLRD1 = data$KLRD1
dd1$TLR4 = data$TLR4
dd1$DCK = data$DCK
dd1$C5 = data$C5
dd1$CCL28 = data$CCL28
dd1$RFX5 = data$RFX5
dd1$EIF2AK2 = data$EIF2AK2
model2 <- glm(y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4 
              + C5 + CCL28 + RFX5 + EIF2AK2 + DCK + pathStage, data = dd1, family = 'binomial')
summary(model2)
## 
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + 
##     CD86 + DDX58 + KLRD1 + TLR4 + C5 + CCL28 + RFX5 + EIF2AK2 + 
##     DCK + pathStage, family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4960  -0.8346   0.2168   0.7899   2.1391  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.73513    0.48444   3.582 0.000341 ***
## MC1R         -0.37369    0.32273  -1.158 0.246901    
## AKT1         -0.41097    0.27016  -1.521 0.128214    
## PAK2          0.22909    0.20232   1.132 0.257502    
## SOS1         -0.26932    0.18330  -1.469 0.141755    
## B2M          -0.40808    0.26533  -1.538 0.124050    
## sex          -0.26646    0.33268  -0.801 0.423165    
## age_c        -0.03155    0.01228  -2.569 0.010209 *  
## CD86          0.37462    0.32247   1.162 0.245341    
## DDX58        -0.21035    0.22115  -0.951 0.341527    
## KLRD1         0.14931    0.25890   0.577 0.564135    
## TLR4          0.30899    0.24808   1.246 0.212944    
## C5            0.19427    0.16974   1.145 0.252397    
## CCL28         0.06700    0.19339   0.346 0.728989    
## RFX5          0.30399    0.20599   1.476 0.140020    
## EIF2AK2       0.32518    0.22064   1.474 0.140542    
## DCK           0.46457    0.22625   2.053 0.040039 *  
## pathStage2   -1.23248    0.56163  -2.194 0.028202 *  
## pathStage3   -2.31474    0.53271  -4.345 1.39e-05 ***
## pathStage4  -17.61213  882.74403  -0.020 0.984082    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 331.89  on 239  degrees of freedom
## Residual deviance: 234.58  on 220  degrees of freedom
## AIC: 274.58
## 
## Number of Fisher Scoring iterations: 13
model2.diag <- glm.diag(model2)
data.frame("cooks" = model2.diag$cook, "leverage" = model2.diag$h) %>%
  pivot_longer(everything(), names_to = "measure") %>%
  mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
  ggplot(aes(x = seq(0.5,240,by=0.5))) +
  geom_line(aes(y = value, col = measure)) +
  geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
  theme_classic() +
  facet_wrap(~measure, scales = "free_y", nrow = 2) +
  scale_color_manual(values = c("orange","cornflowerblue")) +
  labs(x = "observation")

p_res2 <- residuals(model2, type = "pearson")
d_res2 <- residuals(model2, type = "deviance")
data.frame("pearson" = p_res2, "deviance" = d_res2) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation")

data.frame("pearson" = model2.diag$rp, "deviance" = model2.diag$rd) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
  labs(x = "observation", y = "standardized values")

#kable(data.frame("pearson" = p_res2, "deviance" = d_res2,"leverage" = model2.diag$h, "cooks" = model2.diag$cook,"p_standard" = model2.diag$rp, "d_standard" = model2.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
# model 3: baseline + selected top genes + immune response genes + age:immunse response genes
model3 <- glm(y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4
                + age_c:CD86 + age_c:DDX58 + age_c:KLRD1 + age_c:TLR4 , data = dd1, family = 'binomial')
summary(model3)
## 
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + 
##     CD86 + DDX58 + KLRD1 + TLR4 + age_c:CD86 + age_c:DDX58 + 
##     age_c:KLRD1 + age_c:TLR4, family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1112  -0.9507   0.3760   0.9613   2.0786  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.291027   0.202136   1.440  0.14993   
## MC1R        -0.375414   0.295563  -1.270  0.20403   
## AKT1        -0.459284   0.228609  -2.009  0.04453 * 
## PAK2         0.290848   0.176836   1.645  0.10002   
## SOS1        -0.126717   0.159818  -0.793  0.42785   
## B2M         -0.182347   0.220100  -0.828  0.40740   
## sex         -0.317183   0.311342  -1.019  0.30832   
## age_c       -0.029922   0.010791  -2.773  0.00556 **
## CD86         0.295873   0.290758   1.018  0.30887   
## DDX58        0.282801   0.210799   1.342  0.17974   
## KLRD1        0.053540   0.253813   0.211  0.83293   
## TLR4         0.303441   0.233045   1.302  0.19289   
## age_c:CD86  -0.006308   0.020972  -0.301  0.76358   
## age_c:DDX58 -0.047733   0.019392  -2.461  0.01384 * 
## age_c:KLRD1 -0.006409   0.016353  -0.392  0.69510   
## age_c:TLR4   0.022713   0.018525   1.226  0.22017   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 331.89  on 239  degrees of freedom
## Residual deviance: 266.87  on 224  degrees of freedom
## AIC: 298.87
## 
## Number of Fisher Scoring iterations: 5
model3.diag <- glm.diag(model3)
data.frame("cooks" = model3.diag$cook, "leverage" = model3.diag$h) %>%
  pivot_longer(everything(), names_to = "measure") %>%
  mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
  ggplot(aes(x = seq(0.5,240,by=0.5))) +
  geom_line(aes(y = value, col = measure)) +
  geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
  theme_classic() +
  facet_wrap(~measure, scales = "free_y", nrow = 2) +
  scale_color_manual(values = c("orange","cornflowerblue")) +
  labs(x = "observation")

p_res3 <- residuals(model3, type = "pearson")
d_res3 <- residuals(model3, type = "deviance")
data.frame("pearson" = p_res3, "deviance" = d_res3) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation")

data.frame("pearson" = model3.diag$rp, "deviance" = model3.diag$rd) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
  labs(x = "observation", y = "standardized values")

#kable(data.frame("pearson" = p_res3, "deviance" = d_res3,"leverage" = model3.diag$h, "cooks" = model3.diag$cook,"p_standard" = model3.diag$rp, "d_standard" = model3.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
library(VGAM)
## Warning: package 'VGAM' was built under R version 4.1.3
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following objects are masked from 'package:boot':
## 
##     logit, simplex
## The following object is masked from 'package:car':
## 
##     logit
## The following object is masked from 'package:tidyr':
## 
##     fill
# model 5: proportional odds model
dd1$CD86 = data$CD86
dd1$DDX58 = data$DDX58
dd1$KLRD1 = data$KLRD1
dd1$TLR4 = data$TLR4
dd1$pathStage = as.factor(data$stage)
model5 <- vglm(pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4, data = dd1, family = cumulative(parallel=T))
## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()
summary(model5)
## 
## Call:
## vglm(formula = pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + 
##     sex + age_c + CD86 + DDX58 + KLRD1 + TLR4, family = cumulative(parallel = T), 
##     data = dd1)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -1.204631   0.155436  -7.750 9.19e-15 ***
## (Intercept):2  0.288234   0.138098   2.087   0.0369 *  
## (Intercept):3  5.690631   1.008460   5.643 1.67e-08 ***
## MC1R          -0.318350   0.139875  -2.276   0.0228 *  
## AKT1           0.283519   0.124273   2.281   0.0225 *  
## PAK2           0.215181   0.104077   2.068   0.0387 *  
## SOS1           0.059470   0.098668   0.603   0.5467    
## B2M            0.018327   0.127866   0.143   0.8860    
## sex           -0.179287   0.184433  -0.972   0.3310    
## age_c         -0.009620   0.005756  -1.671   0.0946 .  
## CD86          -0.291645   0.162506  -1.795   0.0727 .  
## DDX58          0.245220   0.099962   2.453   0.0142 *  
## KLRD1          0.050054   0.136286   0.367   0.7134    
## TLR4          -0.030667   0.115717  -0.265   0.7910    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), 
## logitlink(P[Y<=3])
## 
## Residual deviance: 504.6351 on 706 degrees of freedom
## 
## Log-likelihood: -252.3176 on 706 degrees of freedom
## 
## Number of Fisher scoring iterations: 14 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):3'
## 
## 
## Exponentiated coefficients:
##      MC1R      AKT1      PAK2      SOS1       B2M       sex     age_c      CD86 
## 0.7273479 1.3277943 1.2400857 1.0612743 1.0184961 0.8358656 0.9904259 0.7470341 
##     DDX58     KLRD1      TLR4 
## 1.2779020 1.0513284 0.9697985
model6 <- vglm(pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4
                + age_c:CD86 + age_c:DDX58 + age_c:KLRD1 + age_c:TLR4 , data = dd1, family = cumulative(parallel = T))
## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()
summary(model6)
## 
## Call:
## vglm(formula = pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + 
##     sex + age_c + CD86 + DDX58 + KLRD1 + TLR4 + age_c:CD86 + 
##     age_c:DDX58 + age_c:KLRD1 + age_c:TLR4, family = cumulative(parallel = T), 
##     data = dd1)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -1.239885   0.163185  -7.598 3.01e-14 ***
## (Intercept):2  0.325951   0.144582   2.254 0.024168 *  
## (Intercept):3  5.822924   1.017363   5.724 1.04e-08 ***
## MC1R          -0.386994   0.142384  -2.718 0.006569 ** 
## AKT1           0.278968   0.128804   2.166 0.030324 *  
## PAK2           0.212728   0.107358   1.981 0.047539 *  
## SOS1           0.069787   0.102485   0.681 0.495906    
## B2M            0.049921   0.137289   0.364 0.716143    
## sex           -0.221438   0.193277  -1.146 0.251918    
## age_c         -0.016624   0.006295  -2.641 0.008274 ** 
## CD86          -0.198693   0.168893  -1.176 0.239420    
## DDX58          0.338315   0.115911   2.919 0.003514 ** 
## KLRD1         -0.058325   0.156900  -0.372 0.710092    
## TLR4          -0.090631   0.125003  -0.725 0.468431    
## age_c:CD86     0.025759   0.013054   1.973 0.048472 *  
## age_c:DDX58   -0.034658   0.010179  -3.405 0.000662 ***
## age_c:KLRD1   -0.027954   0.010426  -2.681 0.007335 ** 
## age_c:TLR4    -0.015682   0.010922  -1.436 0.151045    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), 
## logitlink(P[Y<=3])
## 
## Residual deviance: 489.3675 on 702 degrees of freedom
## 
## Log-likelihood: -244.6837 on 702 degrees of freedom
## 
## Number of Fisher scoring iterations: 15 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):3'
## 
## 
## Exponentiated coefficients:
##        MC1R        AKT1        PAK2        SOS1         B2M         sex 
##   0.6790952   1.3217656   1.2370476   1.0722794   1.0511879   0.8013654 
##       age_c        CD86       DDX58       KLRD1        TLR4  age_c:CD86 
##   0.9835135   0.8198018   1.4025826   0.9433435   0.9133542   1.0260937 
## age_c:DDX58 age_c:KLRD1  age_c:TLR4 
##   0.9659355   0.9724329   0.9844399